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Abstract. It is found experimentally that the coexistence region of a vapor-liquid system or a binary 
mixture is substantially narrowed when the fluid is confined in a aerogel with a high degree of porosity 
(e.g. of the order of 95% to 99%). A Hamiltonian model for this system has recently been introduced 
[1]. We have performed Monte-Carlo simulations for this model to obtain the phase diagram for the 
model. We use a periodic fractal structure constructed by diffusion-limited cluster-cluster aggregation 
(DLCA) method to simulate a realistic gel environment. The phase diagram obtained is qualitatively 
similar to that observed experimentally. We also have observed some metastable branches in the phase 
diagram which have not been seen in experiments yet. These branches, however, might be important in 
the context of recent theoretical predictions and other simulations. 
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1. Introduction 

When a simple liquid or a binary mixture is con- 
fined in a porous material which has a very low 
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density (1-5%) of spatially fixed impurities, such 
as in an aerogel, the coexistence region in the 
phase diagram is substantially narrowed. This 
result has been obtained in a broad class of ex- 
perimental studies, such as vapor-liquid coexis- 
tence of 4 He[2] and Nitrogen[3], binary mixtures 
of isobutyric acid-water[4] and 3 He- 4 He[5], etc. 
In all of these studies, the coexistence curve was 
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shown to change dramatically when the system 
was confined in a low concentration silica aerogel. 

Recent theoretical efforts have been aimed to 
understand the above mentioned behavior. These 
include mean-field type studies of the Random 
Field Ising model[6], a liquid state approach using 
the Ornstein-Zernike equations [7] and numerical 
simulations of a modified version of the Blume- 
Emery-Griffiths model [8]. A very successful ap- 
proach was initiated by Donley and Liu[l]. In 
this reference, the authors introduce a free energy 
functional that takes into account the interactions 
that arise from the contact between the system 
molecules and the aerogel. By minimizing this 
free energy they obtain a coexistence curve which 
is in rough qualitative agreement with the exper- 
imental results. Moreover, the authors go beyond 
this mean field type approach by using a para- 
metric form of the equation of state, combined 
with linear interpolation techniques. Although 
this new approach yields better results than the 
previous mean field treatment, it is not conclusive 
since other parametric models may give different 
results. Moreover, as the authors point out cor- 
rectly, it is very important to study the role of the 
fluctuations. 

In this paper, we go beyond the mean field ap- 
proach and numerically determine the phase di- 
agram of the model introduced in [1] by using 
Monte Carlo methods. In this model, one consid- 
ers a scalar field m(r) and writes down a Hamilto- 
nian which includes bulk terms plus surface terms 
coming from the interaction with the aerogel: 
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-m 2 (r) + ^m 4 (r) 
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The bulk terms, the first volume integral, is 
the usual Ginzburg-Landau model for a scalar 
concentration field m(r) used in binary phase- 
transitions. The additional term given by the sur- 
face integral represents the superficial stress [9] in 
the neighborhood of gel. Here, the volume V is the 
available volume for the fluid and the surface S is 
the set of fluid points in contact with the gel. The 



parameters for this model are: 9, which is related 
to the temperature; \ which sets the width of the 
coexistence curve; the external field (playing the 
role of the chemical potential) H; the surface field 
Hi; and the surface enhancement parameter Q. 

We have performed Monte-Carlo simulations of 
the lattice version of the above Hamiltonian in 
order to find its phase diagram. We consider a 
three-dimensional simple cubic lattice with peri- 
odic boundary conditions. In this lattice, we sim- 
ulate the presence of the aerogel by considering 
that Nq out of the L 3, lattice sites belong to a 
gel structure generated in a way to be explained 
in detail later. We call these sites "gel sites" . In 
the remaining sites (the "field sites" ) we consider 
the scalar variable mi (i = 1, . . . , TV = L 3 — Nq) 
representing the fluid density field. The gradient 
term of eq.(l) is discretized in the usual way: 



|Vm(r)| 2 ^^(m 4 



rm 
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where (ii,i2)*3) stands for the set of right 
nearest-neighbors sites to site i. However, the 
presence of the gel has the effect that in this ex- 
pression for the gradient: only those neighbor sites 
which are actually field sites contribute to the 
sum. Accordingly, we introduce a set of indexes, 
Oi^ defined to be equal to 1, when the site is a 
field site, or 0, when the site is a gel site. The 
gradient term becomes then: 



|Vm(r)| 2 



rm 



(3) 



For the sake of clarity in notation we have ordered 
the N field points such that, from 1 to Nb we 
have "pure bulk" field sites (i.e. those which are 
not in contact with the gel) and from Nb + 1 to 
N we have N$ "surface" field sites, (i.e. those 
in contact with the gel, N s = N - N B ). With 
this convection in mind, the lattice version of the 
Hamiltonian (eq.(l)) can be written as: 



N 



H = ^ [tmf + um\ — hrrii 
i=i 

1 3 
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Fig. 1 . Three dimensional gel structure with a concentra- 
tion of 4%, generated by a realization of the DLCA process, 
in a lattice with L = 32 and periodic boundary conditions. 

N 

+ [-himi+gm?] (4) 

i=N B +l 

Where, t, u, h, hi and g are parameters ob- 
tained by suitable rescaling of the continuum val- 
ues 9, Xi Hi H\, G, respectively. The gel sites in 
this lattice form a periodic fractal structure gen- 
erated by a diffusion - limited-cluster-aggregation 
(DLCA) process[10, 11], which mimics the aggre- 
gation process that form silica gels. The algorithm 
proceeds as follows [12]: 

Let us consider the starting configuration of the 
gel as a collection of aggregates (clusters) contain- 
ing one particle each, the total number of particles 
is No- At a later time, one obtains a collection of 
N a aggregates, the i-th aggregate containing rii 
gel particles, so that 

Y^n i =N G (5) 

i=l 

The aggregates evolve in the following way: an 
aggregate i is chosen at random according to a 
probability p ni which depends on the number of 
particles rn that it contains, given by 

with a = —0.55. Then a space direction is cho- 
sen at random among the six possible directions 
and the cluster is moved by one lattice step in 



that direction (we use periodic boundary condi- 
tions). If the cluster does not collide with any 
other cluster the algorithm continues by choosing 
again another cluster at random and moving it. If 
instead a collision occurs, the two colliding clus- 
ters merge into a new cluster formed by sticking 
together the colliding clusters. The process is re- 
peated until a single cluster remains in the system. 
The resulting fractal dimension of the clusters is 
Dp = 1.9±0.1 which is close to the expected value 
Dp = 1/a w 1.78. In Fig.(l), we can see a pic- 
ture of a fractal gel structure obtained using this 
DLCA process. 

2. The method 

We use the average value of the field as the order 
parameter (M): 

w = (^X>) ( 7 ) 

where, for fixed gel structure, averages are per- 
formed with respect to the distribution e~ n . For 
a given gel structure representative configura- 
tions are obtained by the use of the Monte-Carlo 
method applied to the lattice Hamiltonian (3). We 
have used the simple Metropolis algorithm: a field 
value rrii is proposed to change to a new value 
chosen randomly from a uniform distribution in 
(m, — 5 1 m i + 5) for given S. The new value m£ is 
accepted with a probability given by min[l, c~ A7i ] , 
with AH = H' — H is the change in the Hamilto- 
nian implied by the proposed change. The order 
parameter (M) is computed as an average over dif- 
ferent field configurations. An additional average 
has been performed with respect to 10 different 
gel structures. 

To find the phase diagram, i.e. the dependence 
on the "temperature" t of the order parameter 
(M), we take fixed values for the system param- 
eters u, g and hi, and vary the "temperature" 
t. For each value of the temperature t we com- 
pute the hysteresis loop by using the Monte-Carlo 
method varying the external field h from +h° to 
— h° and vice- versa. We first start at a sufficiently 
high value for ft, (see later) and compute (M)°. 
Next, by keeping the same final configuration for 
the density field, the external field is lowered by an 
amount A to h 1 = h° — A and compute the cor- 
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Fig. 2. Hysteresis loops, (M) versus /i, for two cases: The gel case (o) with a concentration of c = 4% and the no gel case 
(+), both for a parameter value t = —1.26. The two vertical lines show the actual value h where the transition take place, 
in each case, h = —0.18 in the gel case, and h = in the no gel case. Note the presence of little steps in the lower branch 
in the gel hysteresis loop. 



responding value for the order parameter (M) 1 . 
Then the field is changed to ft 2 = ft 1 — A 1 and so 
on until we arrive at —ft . The process is reversed 
by increasing in a similar way the external field to 
reach again +ft°. 

In order to determine accurately the hysteresis 
loop, we do not take a constant value for A 1 but 
we take: 



A°- 



A l - 



^J(h l - ft 1 " 1 ) 2 + ((My - (Mp)V 

(8) 

where a is an additional scale control parameter. 
This means that we control the length along the 
hysteresis curve allowing us to have smooth hys- 
teresis curves. Two typical results for the hys- 
teresis loops are shown in Fig. (2) for the cases of 
no gel and a gel filling 4% of the lattice points. In 
the no-gel case, the hysteresis loop is symmetrical 
around ft = and one can read directly the equi- 
librium values for ±(M) by taking the values at 
h = 0. When the gel is present, we determine the 
equilibrium values for (M) by demanding that the 
Gibbs free energy in the two phases is equal. The 
Gibbs free energy can be obtained by integration 
of the general relation[13]: 



(M) = - 



dG 
dh 



(9) 



By integrating along the upper curve of the hys- 
teresis loop we obtain: 



G {1) (h) = G(h )- / (M)dh (10) 
Jh° 

whereas from the lower part of the hysteresis loop: 



G (2) (ft) = G(-ft°) 



h° 



(M)dh (11) 



The equilibrium values for (M) are read from 
the hysteresis loop at the value of the external 
field ft such that G {1 \h) = G^ 2 \h). In order to 
compute those values for the free energy, accord- 
ing to (10) and (11) we need to know the values 
of G(h°) and G(—h°). For this, we use a suffi- 
ciently large value for ft . For such a large exter- 
nal field, the mean field is a good approximation, 
in such a way that the Gibbs free-energy can be 
replaced by just the internal energy H. So we take 
G(±ft°) w H(±h°). We have taken ft = 10. In 
order to check the validity of mean field for this 
value of ft we have compared the resulting av- 
erage (M)° obtained in the simulation with the 
mean field value obtained from minimizing the 
Hamiltonian H for the same value for ft . Both 
results agree within 1%. 
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Fig. 3. The Gibbs free-energy G{h) versus h obtained by integration, Eqs. (10,11), from the hysteresis curves in Fig.(2), 
in the gel case (o) and the no gel case (+). From these curves we deduce the necessary value for the external field h, where 
the first-order transition takes place. 



In Fig. (3) we plot the results of the numerical 
integrations (10) and (11) both for the no-gel case 
and for one gel configuration of porosity 96% for 
a value of the parameter t — —1.26. As expected, 
in the no gel case, G^\h) and G {2 \h) coincide 
for h = 0. In the gel case, we read from this curve 
the corresponding value for h w —0.18. Using this 
value, we obtain from the upper and lower curves 
of the hysteresis loops, see Fig. (2), the correspond- 
ing values for (M). 

We have found that this method can be used ef- 
ficiently far enough from the critical point. Near 
the critical point, the numerical errors produce 
a large uncertainty in the numerical integrations 
and it is difficult to accurately determine the re- 
quired value of the external field. In those cases, 
we have taken simply an average of the lower and 
upper branches of the hysteresis loop as the values 
for (M). For temperatures above the critical one, 
there is no hysteresis loop. 



3. Phase Diagram 

We present in this section the phase diagram as 
a function of the parameter t for three different 
cases: (i) the no-gel situation, (ii) a gel case with a 



porosity of 88% and (iii) a gel case with a porosity 
of 96%. 

We use in all the cases a lattice with L 3 = 32 3 
sites and the common Hamiltonian parameter val- 
ues u — 0.5, hi = 4 and g = 1. For the factors 
controlling the step size for the variation of the 
external field we take a = 3, A = 0.05. The ini- 
tial values for the hysteresis loop is h° — 10. In 
the gel cases, we have taken averages with respect 
to 10 gel structures. By following the method de- 
scribed in the previous section, we obtain for each 
temperature two values for the order parameter 
(M). These are plotted in Fig. (4) and Fig. (5) for 
a porosity of 88% and 96%, respectively. In these 
figures, we can see clearly the narrowing of the 
coexistence region when the gel is present. For 
smaller porosity (larger fraction of the gel) the 
narrowing is more pronounced as observed in the 
experiments and in accordance with the results 
of mean field theory. Although, as mentioned at 
the end of the previous section, numerical errors 
become large near the critical point, we conclude 
from the figures that the critical temperature is 
lowered when the gel is present. Again, the re- 
duction in the critical temperature is larger for 
smaller porosity. 

However, we note that in the simulations where 
the gel is present, we have found some steps in the 
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Fig. 4- Phase diagrams, in the gel (o), c = 12%, and no gel (continuous line) cases. 



lower curve of the hysteresis loop, as we can see 
for instance in Fig. (2), which could be interpreted 
as signaling a second transition. We have plot- 
ted in the phase diagram additional points corre- 
sponding to the steps in the hysteresis loops. The 
location of those points depends strongly on the 
particular realization of the DLCA process to gen- 



erate the fractal gel structure. This shows up in 
the large error bars for these points in the Fig. (5). 

The phase diagram obtained in this paper is 
qualitatively similar to that observed experimen- 
tally [2]: the coexistence region in presence of gel is 
narrowed and shifted with respect to the non-gel 
situation. There are marks of a second transition 
which also show up in the mean field studies of 
reference ([!]) and also in other simulations of the 
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Fig. 5. Phase diagrams, in the gel (o), c = 4%, and no gel (continuous line) cases. Note in the gel case that we have 
included some points (+), corresponding to the steps found in the hysteresis loops. 
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Lennard-Jones fluid[14], although it has not been 
reported in experimental studies. 
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